function solve_counter = solve_counter(fp,param)
alpha = param(1:11);
TFP = param(12:13) ;  

init_draws      = cell(4,1);
peers_forward   = cell(4,1);
peers_backward  = cell(4,1);
PS_forward      = cell(4,1);

h = 1; %solve for the equilibrium first;

for j = [fp.origin_neigh fp.receiving_neigh] %solve equilibrium only for sending and receiving neighborhoods;

    init_draws{j}     = fp.schools_distrib(j,1) + fp.schools_distrib(j,2).*(fp.init_draws{j}(1:fp.schools_distrib(j,3),1,h) - mean(fp.init_draws{j}(1:fp.schools_distrib(j,3),1,h)))./std(fp.init_draws{j}(1:fp.schools_distrib(j,3),1,h));

    %Set the vector of Shocks for both sending and receiving children;
    if j==fp.origin_neigh

        peers_forward{j}  = fp.peers_forward{j}(1:fp.schools_distrib(j,3),1:fp.schools_distrib(j,3),h,:);
        peers_backward{j} = fp.peers_backward{j}(:,1:fp.schools_distrib(j,3), :,:);
        PS_forward{j}     = fp.PS_forward{j}(1:fp.schools_distrib(j,3),h,:);

    else

        peers_forward{j}  = fp.peers_forward{j}(1:(fp.schools_distrib(j,3)+fp.n_kids_moved),1:(fp.schools_distrib(j,3)+fp.n_kids_moved),h,:);
        peers_backward{j} = fp.peers_backward{j}(:,1:(fp.schools_distrib(j,3)+fp.n_kids_moved), :,:);
        PS_forward{j}     = fp.PS_forward{j}(1:(fp.schools_distrib(j,3)+fp.n_kids_moved),h,:);


    end
end

%Ordering children to be moved from the sending neighborhoods;
temp_init_draws = sort(init_draws{fp.origin_neigh}(:,1), 1, 'ascend');

if fp.do_counter_type==1 %First policy (moving one child);

    if fp.n_kids_moved ==1
        %Moving the Median Child;
        moved_children_skills = temp_init_draws( round( fp.schools_distrib(fp.origin_neigh,3)/2 ) , 1);
    end

    if fp.n_kids_moved==0 %if no one is moved;
        moved_children_skills  = [];
    end
end %end first policy;


if fp.do_counter_type==2 || fp.do_counter_type==3 %Second policy (moving many children);
    %Gradually moving children;
    moved_children_skills = temp_init_draws(fp.shock_child_moved(:,h) , 1 ) ;
end %end second policy;



%%%%%
%Gathering/Ordering Data for Sent and Receiving Children in the Counterfactual Economy;
%%%%%

%Initial Draws (Keep the same shock for every child Baseline vs
%Counterfactual);
fp_parfor.init_draws = [ init_draws{fp.receiving_neigh}; moved_children_skills ];
fp_parfor.peers_forward  = peers_forward{fp.receiving_neigh};
fp_parfor.peers_backward = peers_backward{fp.receiving_neigh};
fp_parfor.PS_forward     = PS_forward{fp.receiving_neigh};

%Indicator if the child is moved or receiving;
fp_parfor.moved_children = [ zeros(length(init_draws{fp.receiving_neigh}),1) ; ones( length(moved_children_skills) , 1)];

fp_parfor.n_simulation = length(fp_parfor.init_draws);

epsilon = 5;
toll= 0.05;

record_epsilon=[3 epsilon];

epsilon_skills = 5;

i = 1;
j = fp.receiving_neigh;


if fp.do_counter_type < 3   %Those are the two policies where we move children and parents change behaviors (fp.do_counter_type=3 we keep parental behavior fixed);

    while (epsilon > toll  && epsilon_skills>toll)  %Fix point;

        fp_parfor.i = i;  %iteration;
        fp_parfor.j = j;  %neighborhood;

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
        %Simulate Forward;
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

        if i == 1

            %first iteration: set the initial guess about the distribution of
            %skills;

            fp_parfor.distrib_PS = zeros( size( fp_parfor.init_draws , 1 ) ,fp.periods-1);

            fp_parfor.inv0{1} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);
            fp_parfor.inv0{2} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);
            fp_parfor.inv0{3} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);

            fp_parfor.inv1{1} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);
            fp_parfor.inv1{2} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);
            fp_parfor.inv1{3} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);

            distrib_skills1 = exp(fp_parfor.init_draws(:,1));
            distrib_skills2 = technology(1,alpha, distrib_skills1 , fp.Max_Time, fp_parfor.distrib_PS(:,1) ,  distrib_skills1 , TFP) ;
            distrib_skills3 = technology(2,alpha, distrib_skills2 , fp.Max_Time, fp_parfor.distrib_PS(:,2) ,  distrib_skills2 , TFP)  ;
            distrib_skills4 = technology(3,alpha, distrib_skills3 , fp.Max_Time, fp_parfor.distrib_PS(:,3) ,  distrib_skills3 , TFP)  ;

            fp_parfor.distrib_skills(:,1) = distrib_skills1;
            fp_parfor.distrib_skills(:,2) = distrib_skills2;
            fp_parfor.distrib_skills(:,3) = distrib_skills3;
            fp_parfor.distrib_skills(:,4) = distrib_skills4;


            fp_parfor.grid{1} = prctile( distrib_skills1 , fp.percentiles_kid ) ;
            fp_parfor.grid_peers{1} = fp_parfor.grid{1} ;
            fp_parfor.grid{2} = prctile( distrib_skills2 , fp.percentiles_kid ) ;
            fp_parfor.grid_peers{2} = fp_parfor.grid{2} ;
            fp_parfor.grid{3} =  prctile(distrib_skills3 , fp.percentiles_kid );
            fp_parfor.grid_peers{3} = fp_parfor.grid{3};

        else
            %solve forward to simulate the distribution of skills and investments;
            sf = solve_forward_counter(param,fp,fp_parfor,h);
            previous_iter_distrib_skills = fp_parfor.distrib_skills;

            fp_parfor.distrib_skills = sf.distrib_skills;
            fp_parfor.distrib_PS = sf.distrib_PS;

            fp_parfor.grid{1} = linspace( prctile(exp(fp_parfor.init_draws(:,1)) , 5) , prctile(exp(fp_parfor.init_draws(:,1)), 95), length(fp.percentiles_kid) );
            fp_parfor.grid_peers{1} = fp_parfor.grid{1} ;
            fp_parfor.grid{2} = linspace( prctile( fp_parfor.distrib_skills(:,2) , 5) , prctile(fp_parfor.distrib_skills(:,2), 95 ), length(fp.percentiles_kid) );
            fp_parfor.grid_peers{2} = fp_parfor.grid{2} ;
            fp_parfor.grid{3} = linspace( prctile( fp_parfor.distrib_skills(:,3), 5 ) , prctile(fp_parfor.distrib_skills(:,3), 95 ), length(fp.percentiles_kid) );
            fp_parfor.grid_peers{3} = fp_parfor.grid{3} ;

            fp_parfor.inv0{1} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);
            fp_parfor.inv0{2} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);
            fp_parfor.inv0{3} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);

            fp_parfor.inv1{1} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);
            fp_parfor.inv1{2} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);
            fp_parfor.inv1{3} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);


        end

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
        %Solve Backward;
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

        sb= solve_backward_counter(param, fp,fp_parfor) ;

        fp_parfor.Hc_tp1_guess_prime = sb.Hc_tp1_guess_prime;

        if i > 1
            eps1=zeros(fp.periods-1,1);
            eps2=zeros(fp.periods-2,1);
            eps_skills=zeros(fp.periods-2,1);


            for t = 1:fp.periods-1

                eps1(t) = max(abs(fp_parfor.estim_policy0{t}(:) - sb.estim_policy0{t}(:)));

                if t==3

                end

                if t < fp.periods-1

                    eps2(t) = max(abs(fp_parfor.estim_policy1{t}(:) - sb.estim_policy1{t}(:)));

                end

                if t > 1
                    deciles_skills_previous_iter = quantile(log(previous_iter_distrib_skills(:,t)),[0.05:0.1:0.95]);
                    deciles_skills_current_iter = quantile(log(sf.distrib_skills(:,t)),[0.05:0.1:0.95]);

                    eps_skills(t) =  max(abs(deciles_skills_previous_iter-deciles_skills_current_iter));

                end


            end

            epsilon = max([eps1;eps2]);
            epsilon_skills = max(eps_skills);
            %keep track of epsilon records;
            ind_rec = mod(i-1,2) + 1;
            record_epsilon(ind_rec) = epsilon;


        end

        fp_parfor.estim_policy0 = sb.estim_policy0;
        fp_parfor.estim_policy1 = sb.estim_policy1;

        fp_parfor.estim_value0 = sb.estim_value0;
        fp_parfor.estim_value1 = sb.estim_value1;


        fprintf('%d ', i)

        i = i + 1;
        %record_epsilon
    end

    fprintf('Solved Neighborhood %d \n ', j)

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Simulate forward once the fixed point is achieved and create the simulated
%dataset;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
fprintf('Simulating Neighborhood %d:\n ', j )

for h = 1 : fp.n_rip_sim %loop over various simulations;

    if h ==fp.n_rip_sim
        fprintf('Simulated Neighborhood %d \n ', j)
    end

    for j = [fp.origin_neigh fp.receiving_neigh] %loop over sending and receiving neighborhoods;

        init_draws{j}     = fp.schools_distrib(j,1) + fp.schools_distrib(j,2).*(fp.init_draws{j}(1:fp.schools_distrib(j,3),1,h) - mean(fp.init_draws{j}(1:fp.schools_distrib(j,3),1,h)))./std(fp.init_draws{j}(1:fp.schools_distrib(j,3),1,h));

        %Shocks;
        if j==fp.origin_neigh

            peers_forward{j}  = fp.peers_forward{j}(1:fp.schools_distrib(j,3),1:fp.schools_distrib(j,3),h,:);
            PS_forward{j}     = fp.PS_forward{j}(1:fp.schools_distrib(j,3),h,:);

        else

            peers_forward{j}  = fp.peers_forward{j}(1:(fp.schools_distrib(j,3)+fp.n_kids_moved),1:(fp.schools_distrib(j,3)+fp.n_kids_moved),h,:);
            PS_forward{j}     = fp.PS_forward{j}(1:(fp.schools_distrib(j,3)+fp.n_kids_moved),h,:);


        end
    end


    temp_init_draws = sort(init_draws{fp.origin_neigh}(:,1), 1, 'ascend');


    if fp.do_counter_type==1

        if fp.n_kids_moved ==1
            %Moving the Median Child;
            moved_children_skills = temp_init_draws( round( fp.schools_distrib(fp.origin_neigh,3)/2 ) , 1);
        end

        if fp.n_kids_moved==0
            moved_children_skills  = [];
        end
    end


    if fp.do_counter_type==2 || fp.do_counter_type==3

        %Gradually moving children;
        moved_children_skills = temp_init_draws(fp.shock_child_moved(:,h) , 1 ) ;


        if fp.do_counter_type==3 %Fixed parental behavior;

            fp_parfor.j = j;

            data_origin_initial = fp.simul_data_baseline{fp.origin_neigh}( fp.simul_data_baseline{fp.origin_neigh}( : ,fp.ind_data.child_age,h)==9 ,:,h);
            data_origin_10      = fp.simul_data_baseline{fp.origin_neigh}( fp.simul_data_baseline{fp.origin_neigh}( : ,fp.ind_data.child_age,h)==10 ,:,h);
            data_origin_11      = fp.simul_data_baseline{fp.origin_neigh}( fp.simul_data_baseline{fp.origin_neigh}( : ,fp.ind_data.child_age,h)==11 ,:,h);
            temp_data           = sortrows( [data_origin_initial(:,[fp.ind_data.child_C   fp.ind_data.inv    fp.ind_data.ParStyle])  data_origin_10(:,[ fp.ind_data.inv    fp.ind_data.ParStyle]) data_origin_11(:,[ fp.ind_data.inv    fp.ind_data.ParStyle])   ] , 1) ;

            data_receiving_inv  =[  fp.simul_data_baseline{fp.receiving_neigh}( fp.simul_data_baseline{fp.receiving_neigh}( : ,fp.ind_data.child_age,h)==9  ,fp.ind_data.inv,h), ...
                fp.simul_data_baseline{fp.receiving_neigh}( fp.simul_data_baseline{fp.receiving_neigh}( : ,fp.ind_data.child_age,h)==10 ,fp.ind_data.inv,h), ...
                fp.simul_data_baseline{fp.receiving_neigh}( fp.simul_data_baseline{fp.receiving_neigh}( : ,fp.ind_data.child_age,h)==11 ,fp.ind_data.inv,h) ];

            data_receiving_PS   =[  fp.simul_data_baseline{fp.receiving_neigh}( fp.simul_data_baseline{fp.receiving_neigh}( : ,fp.ind_data.child_age,h)==9  ,fp.ind_data.ParStyle,h), ...
                fp.simul_data_baseline{fp.receiving_neigh}( fp.simul_data_baseline{fp.receiving_neigh}( : ,fp.ind_data.child_age,h)==10 ,fp.ind_data.ParStyle,h), ...
                fp.simul_data_baseline{fp.receiving_neigh}( fp.simul_data_baseline{fp.receiving_neigh}( : ,fp.ind_data.child_age,h)==11 ,fp.ind_data.ParStyle,h) ];

            fp_parfor.fixed_inv = [ data_receiving_inv ; temp_data( fp.shock_child_moved(:,h) , [ 2 4 6 ] )];
            fp_parfor.fixed_PS  = [ data_receiving_PS  ; temp_data( fp.shock_child_moved(:,h) , [ 3 5 7 ] )];

        end


    end
    %collect children who have been moved;
    id_moved_children(:,h) = ismember( init_draws{fp.origin_neigh}(:,1) , moved_children_skills);


    %%%%%
    %Gathering/Ordering Data for Sent and Receiving Children in the Counterfactual Economy;
    %%%%%

    %Initial Draws (Keep the same shock for every child Baseline vs
    %Counterfactual);
    fp_parfor.init_draws = [ init_draws{fp.receiving_neigh}; moved_children_skills ];
    fp_parfor.peers_forward  = peers_forward{fp.receiving_neigh};
    fp_parfor.peers_backward = peers_backward{fp.receiving_neigh};
    fp_parfor.PS_forward     = PS_forward{fp.receiving_neigh};

    %Indicator if the child is moved or receiving;
    fp_parfor.moved_children = [ zeros(length(init_draws{fp.receiving_neigh}),1) ; ones( length(moved_children_skills) , 1)];

    fp_parfor.n_simulation = length(fp_parfor.init_draws);

    %simulate economy;
    sf = solve_forward_counter(param,fp,fp_parfor,h);
    %collect simulated data;
    sim_data_counter(:,:,h) = sf.sim_data ;
    sim_data_network_counter(:,:,h) = sf.sim_data_network;



end


solve_counter.sim_data = sim_data_counter;
solve_counter.sim_data_network = sim_data_network_counter;
solve_counter.id_moved_children = id_moved_children;
